---
title: "MEA firing data"
author: "Adam Bindas"
date: "3/24/2022"
output:
  html_document:
   code_folding: hide
---

```{r setup, include=FALSE}
# Load Packages
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
library(gt)
library(tidyverse)
library(ggplot2)
library(pscl)
library(boot)
library(moments)
library(lmtest)
library(dplyr) 
library(ggpubr)
library(GLMMadaptive)
library(glmmTMB)
library(MASS)
library(RColorBrewer)
library(bbmle)
library(emmeans)
library(DHARMa)
library(gap)

# Import Datasets
mea0 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/MEA data analysis/Raw Firing Data 22-06-05.csv", header=TRUE, stringsAsFactors=FALSE)
```

```{r}
# Define Variables
mea.pre <- subset(mea0, ground == 0 & unit.no == 0)
#hist(mea.pre$spike.per.sec, breaks = 2000, xlim=c(0,0.1))
mea.pre$Groups.f <- factor(mea.pre$group, levels = c("control", "pff", "but",  "pff.4x", "pff.4x.but"))
mea.pre$well.no <- factor(mea.pre$well)
mea.pre$cond <- factor(mea.pre$cond, levels = c("spont", "dop","ach", "brainphys", "brainphys+ach","neurobasal","cap","estim","estim+dop"))
mea.pre$day.no <- factor(mea.pre$day)
mea.pre$rat.no <- factor(mea.pre$replicate)
mea.pre$dosage <- factor(mea.pre$dosage.um)
mea.pre$simple <- ifelse(mea.pre$Groups.f == "control" | mea.pre$Groups.f == "pff" & mea.pre$Groups.f != "pff.4x", 1, 0)
mea.pre$high.activity <- ifelse(mea.pre$spike.per.sec > 0.5, 1, 0)
mea.pre$seed.count <- factor(mea.pre$seed.count)
mea.pre$arac <- factor(mea.pre$arac)
mea.pre$postpff <- ifelse(mea.pre$day > 7, 1, 0)
mea.pre$elect.no <- factor(mea.pre$electrode)
mea.pre$w.id <- factor(paste(mea.pre$mea.id,"-",mea.pre$well.no))
mea.pre$unique.id <- factor(paste(mea.pre$mea.id,"-",mea.pre$well.no,"-",mea.pre$electrode))

# Subset Data
mea.glmmtmb <- subset(mea.pre, 
                      ground == 0 & 
                      unit.no == 0 & 
                      arac == 0)

mea.glmmtmb.count <- subset(mea.glmmtmb,
                            (cond == "spont" | cond == "ach" | cond == "dop" | 
                               cond == "estim" | cond == "estim +dop") &
                            (Groups.f == "control" | Groups.f == "pff" | Groups.f == "pff.4x"))

#Determine Sample Size
n.mea <- aggregate(data = mea.glmmtmb.count,
                          rat.no ~ method.no + Groups.f + day.no + cond,
                          function(rat.no) length(unique(rat.no)))
m.mea <- aggregate(data = mea.glmmtmb.count,
                          mea.id ~ method.no + Groups.f + day.no + cond,
                          function(mea.id) length(unique(mea.id)))
w.mea <- aggregate(data = mea.glmmtmb.count,
                          w.id ~ method.no + Groups.f + day.no + cond,
                          function(w.id) length(unique(w.id)))
i.mea <- aggregate(data = mea.glmmtmb.count,
                          unique.id ~ method.no + Groups.f + day.no + cond,
                          function(unique.id) length(unique(unique.id)))

# Merge Sample Size to a table
mea.replicates <- cbind(n.mea,m.mea[,5],w.mea[,5],i.mea[,5])
colnames(mea.replicates) <- c("Method","Group","Day","Condition","N","mea","wells","no.electrodes")
mea.replicates

         
```

```{r glmm inspection}

#hist(mea.glmmtmb$spike.per.sec, breaks = 2000, xlim=c(0,0.1))
# Define variables for quantification
mea.glmmtmb$nspikec <-round(mea.glmmtmb$spike.per.sec*5*60)
mea.glmmtmb$Groups.f <- factor(mea.glmmtmb$group, levels = c("control", "pff", "but",  "pff.4x", "pff.4x.but"))
mea.glmmtmb$well.no <- factor(mea.glmmtmb$well)
mea.glmmtmb$cond <- factor(mea.glmmtmb$cond, levels = c("spont", "dop","ach", "brainphys", "brainphys+ach","neurobasal","cap","estim","estim+dop"))
mea.glmmtmb$day.no <- factor(mea.glmmtmb$day)
mea.glmmtmb$rat.no <- factor(mea.glmmtmb$replicate)
mea.glmmtmb$dosage <- factor(mea.glmmtmb$dosage.um)
mea.glmmtmb$simple <- ifelse(mea.glmmtmb$Groups.f == "control" | mea.glmmtmb$Groups.f == "pff" & mea.glmmtmb$Groups.f != "pff.4x", 1, 0)
mea.glmmtmb$high.activity <- ifelse(mea.glmmtmb$spike.per.sec > 0.025, 1, 0)
mea.glmmtmb$seed.count <- factor(mea.glmmtmb$seed.count)
mea.glmmtmb$arac <- factor(mea.glmmtmb$arac)
mea.glmmtmb$postpff <- ifelse(mea.glmmtmb$day > 7, 1, 0)

# Define Day Names for Graphs
day.names.mea <- c(`10` = "Day 10",
                   `14` = "Day 14", 
                   `15` = "Day 15", 
                   `20` = "Day 20", 
                   `21` = "Day 21", 
                   `25` = "Day 25",
                   `28` = "Day 28")

day.names.mea2 <- c(`5` = "Day 5",
                   `7` = "Day 7", 
                   `10` = "Day 10",
                   `14` = "Day 14", 
                   `15` = "Day 15", 
                   `20` = "Day 20", 
                   `21` = "Day 21", 
                   `25` = "Day 25",
                   `28` = "Day 28")

# Additional subsetting of datasets
mea.glmm.v1 <- subset(mea.glmmtmb,
                     postpff == 1 &
                     cond == "spont" &
                     arac == 0 &
                     day.no != 29 &
                     (Groups.f == "control" | Groups.f == "pff" | Groups.f == "pff.4x"))

mea.glmm.v2 <- subset(mea.glmmtmb, 
                     postpff == 1 &
                     method.no == 1 &
                     (cond == "spont" | cond == "ach" | cond == "dop") &
                     arac == 0 &
                     day.no != 29 &
                     (Groups.f == "control" | Groups.f == "pff" | Groups.f == "pff.4x"))

mea.glmm.v3 <- subset(mea.glmmtmb, 
                     postpff == 1 &
                     method.no == 2 &
                     arac == 0 &
                     day.no != 29 &
                     (Groups.f == "control" | Groups.f == "pff" | Groups.f == "pff.4x"))

# max.spike <- max(mea.glmm1$nspikec)
# mea.glmm1$beta.spike <- mea.glmm1$nspikec/max.spike
# mea.glmm1$beta.adj.spike <- mea.glmm1$beta.spike * 0.999
```

```{r spontaneous histogram}
# Develop histogram of spontaneous collected data
plot.hist2 <- ggplot(data = mea.glmm.v1, aes(x=nspikec)) +
  geom_histogram(alpha=0.2, position="identity", binwidth = 1) +
  ylab("Count") + 
  xlab("Spike Count (Normalized to 5 minute Recording)") + 
  labs(
    title = "Distribution of Spike Count (Overall Spontaneous Activity)") + 
  theme(text = element_text(size = 14)) 

# Crop histogram
plot.hist2 +
  coord_cartesian(xlim=c(0,20)) 

# Subset Dataset for specific histogram
mea.hist.rep <- subset(mea0, ground == 0 & unit.no == 0 &
                     arac == 0 &
                     day == 15 &
                     cond == "spont" &
                     replicate == 21 &
                     (group == "control" | group == "pff" | group == "pff.4x"))
mea.hist.rep$Groups.f <- factor(mea.hist.rep$group)

# Plot histogram for specific timepoint
plot.histrep <- ggplot(data = mea.hist.rep, aes(x=spike.per.sec, fill=Groups.f)) +
  geom_histogram(alpha=0.2, position="identity", binwidth = 0.005) +
  ylab("Count") + 
  xlab("Spikes per second") + 
  labs(
    title = "Representative histogram of spontaneous spike activity",
    subtile = "Culture Day 20",
    fill = "Dosage") +
  scale_fill_discrete(labels=c("Control", "1 µg PFF","4 µg PFF")) +
  theme(text=element_text(size=14))

# Crop histogram
plot.histrep +
  coord_cartesian(xlim=c(0,0.03)) 

# Summarize values
summary(mea.glmm.v1$nspikec > 20)
summary(mea.glmm.v1$nspikec == 0)
```

```{r Method 1 analysis}
#Zero Inflated Models for Model 1 Analysis (Neurotransmitters)
fit.meth1.zi.pois <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family = poisson)

fit.meth1.zi.n.binom.1 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.zi.n.binom.2 <- update(fit.meth1.zi.pois,family=nbinom1)
fit.meth1.zi.n.binom.t1 <- update(fit.meth1.zi.pois,family=truncated_nbinom2)
fit.meth1.zi.n.binom.t2 <- update(fit.meth1.zi.pois,family=truncated_nbinom1)

#No zero inflation
fit.meth1.pois <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v2,
                           ziformula = ~0,
                           family = poisson)

fit.meth1.n.binom.1 <- update(fit.meth1.pois,family=nbinom2)
fit.meth1.n.binom.2 <- update(fit.meth1.pois,family=nbinom1)

#Model Comparison
anova(fit.meth1.zi.pois, 
      fit.meth1.zi.n.binom.1,fit.meth1.zi.n.binom.2
      ,fit.meth1.zi.n.binom.t1,fit.meth1.zi.n.binom.t2,
      fit.meth1.pois, 
      fit.meth1.n.binom.1,fit.meth1.n.binom.2)

anova(fit.meth1.zi.n.binom.1,fit.meth1.zi.n.binom.t1,fit.meth1.zi.n.binom.t2,
      fit.meth1.n.binom.1)

# Identify optimal model
AICtab(fit.meth1.zi.pois, 
      fit.meth1.zi.n.binom.1,fit.meth1.zi.n.binom.2,
      fit.meth1.zi.n.binom.t1,fit.meth1.zi.n.binom.t2,
      fit.meth1.pois, 
      fit.meth1.n.binom.1,fit.meth1.n.binom.2)

# Summarize model fit
summary(fit.meth1.zi.n.binom.1)

# Develop multiple models from optimal distribution
fit.meth1.v1 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v2 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                    (1|rat.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v3 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                    (1|rat.no) + (1|rat.no:mea.id),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v4 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                    (1|elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v5 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v6 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|rat.no:elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v7 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1 + day|elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v8 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1 + day|rat.no) + (1 + day|rat.no:elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v9 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|rat.no:well.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v10 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|rat.no:well) + (1|elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v11 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|rat.no:well.no) + (1|rat.no:well.no:elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth1.v12 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day + seed.count +
                                   (1|rat.no) + (1|rat.no:elect.no),
                           data = mea.glmm.v2,
                           ziformula = ~1,
                           family=nbinom2)

# Idenfity optimal model
anova(fit.meth1.v1,fit.meth1.v2,fit.meth1.v3,fit.meth1.v4,fit.meth1.v5,fit.meth1.v6,fit.meth1.v7,
      fit.meth1.v8,fit.meth1.v9,fit.meth1.v10,fit.meth1.v11,fit.meth1.v12)
anova(fit.meth1.v2,fit.meth1.v4,fit.meth1.v5,fit.meth1.v6,fit.meth1.v10,fit.meth1.v11,fit.meth1.v12)
anova(fit.meth1.v4,fit.meth1.v5,fit.meth1.v6,fit.meth1.v11,fit.meth1.v12)
anova(fit.meth1.v5,fit.meth1.v6,fit.meth1.v12)

# Plot Model Fit
simulationOutput.meth1 <- simulateResiduals(fittedModel = fit.meth1.v6)
plot(simulationOutput.meth1)
```

``` {r method 1 post hoc}
# Emmeans statistical analysis
em.fit.meth1 <- emmeans(fit.meth1.v6, ~Groups.f * cond)
pwpm(em.fit.meth1)
plot(em.fit.meth1)
summary(em.fit.meth1)

#Plots
cond.names.meth1.a <- c(`spont` = "Spontaneous",
                   `dop` = "Dopamine", 
                   `ach` = "Acetylcholine")

cond.names.meth1.b <- c(`control` = "Control",
                   `pff` = "1 µg PFF", 
                   `pff.4x` = "4 µg PFF")

# hist(mea.glmm.v2$nspikec, xlim = c(0,10), breaks = 1000)
# 
# ggplot(data = mea.glmm.v2, aes(x = Groups.f, y = nspikec, fill = cond)) +
#     geom_boxplot() + coord_cartesian(ylim=c(0,50))

plot.em.fit.meth1 <- plot(em.fit.meth1, comparisons = TRUE, plotIT=FALSE) +
  labs(title = "Stimulated Spike Count by Group (Method 1)",
       x = "Spike Count (Normalized to 5-minute Recording)",
       y = "Groups") +
    scale_y_discrete(labels = c(
      "Control Spontaneous", "1 µg PFF Spontaneous","4 µg PFF Spontaneous",
      "Control Dopamine", "1 µg PFF Dopamine","4 µg PFF Dopamine", 
      "Control Acetylcholine", "1 µg PFF Acetylcholine","4 µg PFF Acetylcholine"
    )) + 
  theme(text = element_text(size = 14), axis.title.y = element_blank()) 

plot.em.fit.meth1



#Sample Size determination
n.mea.m1 <- aggregate(data = mea.glmm.v2,
                          rat.no ~ Groups.f + cond ,
                          function(rat.no) length(unique(rat.no)))
m.mea.m1 <- aggregate(data = mea.glmm.v2,
                          mea.id ~ Groups.f + cond ,
                          function(mea.id) length(unique(mea.id)))
w.mea.m1 <- aggregate(data = mea.glmm.v2,
                          w.id ~ Groups.f + cond ,
                          function(w.id) length(unique(w.id)))
i.mea.m1 <- aggregate(data = mea.glmm.v2,
                          unique.id ~ Groups.f + cond ,
                          function(unique.id) length(unique(unique.id)))

# Merge sample size values to a single table
mea.replicates.meth1 <- cbind(n.mea.m1,m.mea.m1[,3],w.mea.m1[,3],i.mea.m1[,3])
colnames(mea.replicates.meth1) <- c("Group","Condition","N","mea","wells","no.electrodes")
mea.replicates.meth1

mea.sample.meth1 <- rbind(
  range(mea.replicates.meth1$N),
  range(mea.replicates.meth1$mea),
  range(mea.replicates.meth1$wells),
  range(mea.replicates.meth1$no.electrodes))
rownames(mea.sample.meth1) <- c("N","mea","wells","no.electrodes")
colnames(mea.sample.meth1) <- c("Minimum","Maximum")
mea.sample.meth1
```

```{r Method 1 voilin plot}
# Develop violin plot of values
violin.meth1.1 <- ggplot(mea.glmm.v2, aes(x=cond, y=nspikec, fill = Groups.f)) + 
  geom_violin() +
  labs(
    title = "Normalized Spike Activity (Method 1)",
    y = "Spike Count (Normalized to 5-minute recording)",
    x = "Stimulation",
    fill = "Group", 
  ) +
  scale_x_discrete(labels=c("Spontaneous", "Dopamine", "Acetylcholine")) +
  theme(text = element_text(size = 14)) +
  scale_fill_viridis_d(labels=c("Control", "1 µg PFF", "4 µg PFF")) +
  geom_signif(
    y_position = c(50,100,150), xmin = c(0.7, 0.7, 1), xmax = c(1.7, 2.7, 3),
    annotation = c("*", "*","***"), tip_length = 0.001, size = 0.8, textsize=8, vjust=0.4
  )

violin.meth1.1

violin.meth1.2 <- ggplot(mea.glmm.v2, aes(x=cond, y=nspikec, fill = Groups.f)) + 
  geom_violin() +
  labs(
    title = "Normalized Spike Activity (Method 1)",
    y = "Spike Count (Normalized to 5-minute recording)",
    x = "Stimulation",
    fill = "Group", 
  ) +
  coord_cartesian(ylim = c(0,20)) +
  scale_x_discrete(labels=c("Spontaneous", "Dopamine", "Acetylcholine")) +
  theme(text = element_text(size = 14)) +
  scale_fill_viridis_d(labels=c("Control", "1 µg PFF", "4 µg PFF")) +
  geom_signif(
    y_position = c(12,14,16), xmin = c(0.7, 0.7, 1), xmax = c(1.7, 2.7, 3),
    annotation = c("*", "*","***"), tip_length = 0.001, size = 0.8, textsize=8, vjust=0.4
  )

violin.meth1.2


```


```{r Method 2 analysis}
#Zero Inflated Distribution Models
fit.meth2.zi.pois <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family = poisson)

fit.meth2.zi.n.binom.1 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.zi.n.binom.2 <- update(fit.meth2.zi.pois,family=nbinom1)
fit.meth2.zi.n.binom.t1 <- update(fit.meth2.zi.pois,family=truncated_nbinom2)
fit.meth2.zi.n.binom.t2 <- update(fit.meth2.zi.pois,family=truncated_nbinom1)

#No zero inflation distribtuion models
fit.meth2.pois <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v3,
                           ziformula = ~0,
                           family = poisson)

fit.meth2.n.binom.1 <- update(fit.meth2.pois,family=nbinom2)
fit.meth2.n.binom.2 <- update(fit.meth2.pois,family=nbinom1)

#Model Comparison to iddentify optimal distribution for dataset
anova(fit.meth2.zi.pois, 
      fit.meth2.zi.n.binom.1,fit.meth2.zi.n.binom.2
      ,fit.meth2.zi.n.binom.t1,fit.meth2.zi.n.binom.t2,
      fit.meth2.pois, 
      fit.meth2.n.binom.1,fit.meth2.n.binom.2)

anova(fit.meth2.zi.n.binom.1,fit.meth2.zi.n.binom.t1,fit.meth2.zi.n.binom.t2,
      fit.meth2.n.binom.1)

AICtab(fit.meth2.zi.pois, 
      fit.meth2.zi.n.binom.1,fit.meth2.zi.n.binom.2,
      fit.meth2.zi.n.binom.t1,fit.meth2.zi.n.binom.t2,
      fit.meth2.pois, 
      fit.meth2.n.binom.1,fit.meth2.n.binom.2)

summary(fit.meth2.zi.n.binom.1)

# Develop multiple models from optimal distribution
fit.meth2.v1 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day,
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v2 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                    (1|rat.no),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v3 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                    (1|rat.no) + (1|rat.no:mea.id),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v4 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                    (1|elect.no),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v5 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|elect.no),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v6 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|rat.no:elect.no),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v7 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1 + day|elect.no),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v8 <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|rat.no:mea.id) + (1|rat.no:mea.id:elect.no),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

fit.meth2.v9 <- glmmTMB(nspikec ~ Groups.f * cond + seed.count +
                                   (1|rat.no) + (1|rat.no:elect.no),
                           data = mea.glmm.v3,
                           ziformula = ~1,
                           family=nbinom2)

# Identify optimal model
anova(fit.meth2.v1,fit.meth2.v2,fit.meth2.v3,fit.meth2.v4,fit.meth2.v5,fit.meth2.v6,fit.meth2.v7,
       fit.meth2.v8,fit.meth2.v9)
anova(fit.meth2.v2,fit.meth2.v4,fit.meth2.v5,fit.meth2.v6,
       fit.meth2.v8,fit.meth2.v9)
anova(fit.meth2.v4,fit.meth2.v5,fit.meth2.v6,
       fit.meth2.v8,fit.meth2.v9)
anova(fit.meth2.v5,fit.meth2.v6,
       fit.meth2.v8,fit.meth2.v9)
anova(fit.meth2.v6,fit.meth2.v8,fit.meth2.v9)
anova(fit.meth2.v8,fit.meth2.v9)

# Plot Model Fit (10,000 simulations to allow for removal of outliers)
simulationOutput.meth2 <- simulateResiduals(fittedModel = fit.meth2.v8, n = 10000)
plot(simulationOutput.meth2)

testOutliers(simulationOutput.meth2, type = "bootstrap")
# length(outliers(simulationOutput.meth2))

# Remove outliers
mea.glmm.v3.2 <- mea.glmm.v3[-c(outliers(simulationOutput.meth2)),]

# Rerun model without outliers
fit.meth2.v8.redo <- glmmTMB(nspikec ~ Groups.f + cond + Groups.f * cond + day +
                                   (1|rat.no) + (1|rat.no:mea.id) + (1|rat.no:mea.id:elect.no),
                           data = mea.glmm.v3.2,
                           ziformula = ~1,
                           family=nbinom2)

# Plot model fit
simulationOutput.meth2.v2 <- simulateResiduals(fittedModel = fit.meth2.v8.redo, n = 10000)
plot(simulationOutput.meth2.v2)
```

``` {r method 2 post analysis}
# Statistical summary
summary(fit.meth2.v8.redo)

em.fit.meth2 <- emmeans(fit.meth2.v8.redo, ~Groups.f * cond)
pwpm(em.fit.meth2)
summary(em.fit.meth2)

#Plots
cond.names.meth2.a <- c(`spont` = "Spontaneous",
                   `estim` = "E-Stim", 
                   `estim+dop` = "E-Stim + Dopamine")

cond.names.meth2.b <- c(`control` = "Control",
                   `pff` = "1 µg PFF", 
                   `pff.4x` = "4 µg PFF")

plot.em.fit.meth2 <- plot(em.fit.meth2, plotIT=FALSE, comparisons = TRUE) +
  labs(title = "Stimulated Spike Count by Group (Method 2)",
       x = "Spike Count (Normalized to 5-minute Recording)",
       y = "Groups") +
    scale_y_discrete(labels = c(
      "Control Spontaneous","1 µg PFF Spontaneous","4 µg PFF Spontaneous",
      "Control E-Stim","1 µg PFF E-Stim","4 µg PFF E-Stim",
      "Control E-Stim + Dopamine","1 µg PFF E-Stim + Dopamine","4 µg PFF E-Stim + Dopamine")) + 
  theme(text = element_text(size = 14), axis.title.y = element_blank()) 

plot.em.fit.meth2

#Sample Size
n.mea.m2 <- aggregate(data = mea.glmm.v3,
                          rat.no ~ Groups.f + cond + day.no,
                          function(rat.no) length(unique(rat.no)))
m.mea.m2 <- aggregate(data = mea.glmm.v3,
                          mea.id ~ Groups.f + cond + day.no,
                          function(mea.id) length(unique(mea.id)))
w.mea.m2 <- aggregate(data = mea.glmm.v3,
                          w.id ~ Groups.f + cond + day.no,
                          function(w.id) length(unique(w.id)))
i.mea.m2 <- aggregate(data = mea.glmm.v3,
                          unique.id ~ Groups.f + cond + day.no,
                          function(unique.id) length(unique(unique.id)))

# Merge sample size into a single table
mea.replicates.meth2 <- cbind(n.mea.m2,m.mea.m2[,4],w.mea.m2[,4],i.mea.m2[,4])
colnames(mea.replicates.meth2) <- c("Group","Condition","Day","N","mea","wells","no.electrodes")
mea.replicates.meth2

mea.sample.meth2 <- rbind(
  range(mea.replicates.meth2$N),
  range(mea.replicates.meth2$mea),
  range(mea.replicates.meth2$wells),
  range(mea.replicates.meth2$no.electrodes))
rownames(mea.sample.meth2) <- c("N","mea","wells","no.electrodes")
colnames(mea.sample.meth2) <- c("Minimum","Maximum")
mea.sample.meth2

violin.meth2.1 <- ggplot(mea.glmm.v2, aes(x=cond, y=nspikec, fill = Groups.f)) + 
  geom_violin() +
  labs(
    title = "Normalized Spike Activity (Method 1)",
    y = "Spike Count (Normalized to 5-minute recording)",
    x = "Stimulation",
    fill = "Group", 
  ) +
  scale_x_discrete(labels=c("Spontaneous", "Dopamine", "Acetylcholine")) +
  theme(text = element_text(size = 14)) +
  scale_fill_viridis_d(labels=c("Control", "1 µg PFF", "4 µg PFF")) +
  geom_signif(
    y_position = c(50,100,150), xmin = c(0.7, 0.7, 1), xmax = c(1.7, 2.7, 3),
    annotation = c("*", "*","***"), tip_length = 0.001, size = 0.8, textsize=8, vjust=0.4
  )

violin.meth2.1

violin.meth2.2 <- ggplot(mea.glmm.v2, aes(x=cond, y=nspikec, fill = Groups.f)) + 
  geom_violin() +
  labs(
    title = "Normalized Spike Activity (Method 1)",
    y = "Spike Count (Normalized to 5-minute recording)",
    x = "Stimulation",
    fill = "Group", 
  ) +
  coord_cartesian(ylim = c(0,20)) +
  scale_x_discrete(labels=c("Spontaneous", "Dopamine", "Acetylcholine")) +
  theme(text = element_text(size = 14)) +
  scale_fill_viridis_d(labels=c("Control", "1 µg PFF", "4 µg PFF")) +
  geom_signif(
    y_position = c(12,14,16), xmin = c(0.7, 0.7, 1), xmax = c(1.7, 2.7, 3),
    annotation = c("*", "*","***"), tip_length = 0.001, size = 0.8, textsize=8, vjust=0.4
  )

violin.meth2.2
```
